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Dear  Dr.  Rood: 


The  purpose  of  this  letter  is  to  transmit  the  eleventh  quarterly  report  for  ONR  Grant  N00014-91-J- 
1271,  "An  Experimental  Study  of  Plunging  Liquid  Jet  Induced  Air  Carryunder  and  Dispersion," 
(Lahey  &  Drew  -  Co-PI). 


During  this  report  period  the  analysis  of  the  pool  surface  depression  by  the  plunging  liquid  jet 
(which  leads  to  the  entrainment  of  air)  was  completed.  This  has  resulted  in  a  technical  paper 
entitled,  "An  Analysis  of  Pool  Surface  Deformation  Due  to  a  Plunging  Liquid  Jet.  This  paper  will 
be  submitted  to  the  Journal  of  Fluid  Mechanics  (JEM)  for  archival  publication. 


We  are  continuing  our  CFD  analysis  of  the  plunging  liquid  jet  data  we  have  taken.  In  particular, 
a  more  generalized  two-fluid  model  is  being  coped.  Moreover,  in  the  future  we  plan  to  analyze  the 
two-phase  shear  flow  data  of  Professor  J.  Lasheras  (UC-SD)  once  it  is  available  to  us. 

If  you  need  any  further  information  concerning  this  report  please  don’t  hesitate  to  contact  me  [(518) 
276-8579]  or  Professor  Drew  [(518)  276-69031. 


Sincerely  yours, 

i)  4  -/ ,  .  /'! 

Dr.  R.T.  Lahey,  jJ. 

The  Edward  E.  Hood,  Jr.  Professor  of  Engineering 
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AN  ANALYSIS  OF  POOL  SURFACE  DEFORMATION 
DUE  TO  A  PLUNGING  LIQUID  JET 


F.  Bonetto 
D.A.  Drew 
R.T.  Lahey,  Jr. 

Center  for  Multiphase  Research 
Rensselaer  Polytechnic  Institute 
Troy,  NY  12180-3590  USA 


ABSTRACT 

When  a  liquid  jet  impacts  a  pool  containing  the  same  liquid  and  surrounded 
by  a  still  gas,  a  surface  depression  is  produced.  The  surface  shape  is  determined  by 

1  2 

2  g  Xo 

the  Weber  number.  We  =  - — —  and  the  Bond  number,  B^^  =  — — —  .  In  this 

work  the  shape  of  the  surface  is  obtained  as  a  function  of  the  Weber  number  and 
Bond  number  by  using  a  non-singular  perturbation  technique. 

INTRODUCTION 

The  entrainment  of  non-condensible  gases  by  a  plunging  liquid  jet  impacting 
a  liquid  pool  is  important  for  some  practical  problems.  For  example,  the  absorption 
of  greenhouse  gases  into  the  ocean  has  been  hypothesized  to  be  highly  dependent 
upon  the  air  carryunder  that  occurs  during  breaking  waves,  which  represent  a  type 
of  plunging  liquid  jet  [Monahan,  1991;  Kerman,  1984].  Other  applications  include 
some  type  of  liquid /gas  chemical  reactors.  In  order  to  enhance  the  reaction  rate,  a  jet 
of  liquid  entrains  the  surrounding  gas  reactant  forming  a  bubbly  two-phase  jet.  The 
reaction  rate  is  enhanced  because  of  the  increase  of  the  interfacial  area  density  and 
the  pseudo-turbulence  produced  by  the  entrained  gas  bubbles. 

Other  less  obvious  applications  are  associated  with  the  dynamics  of  the  slug 
flow  regime  in  two-phase  flow.  For  example,  in  vertical  slug  flow,  the  Taylor 
bubbles  move  upward  with  an  almost  constant  velocity.  The  liquid  film 


QUALITY  INSPECTEI'  8 


2 


surrounding  the  Taylor  bubbles  drains  downward  forming  an  annular  liquid  jet. 
Therefore,  a  gas  entrainment  problem,  similar  to  the  one  studied  in  this  work, 
occurs  behind  the  Taylor  bubbles.  In  particular,  a  Kelvin-Helmholtz  instability 
occurs  at  the  interface  in  the  rear  portion  of  the  Taylor  bubble  such  that  small  gas 
bubbles  are  entrained  into  the  liquid  plug  behind  the  Taylor  bubble. 

DISCUSSION 

Most  prior  theoretical  studies  were  done  for  liquid  jets  with  very  low  liquid 
velocities,  and  had  applications  for  fiber  coating,  etc.  Lezzi  &  Prosperetti  [1991] 
performed  a  stability  analysis  of  the  gas  entrained  by  a  plunging  liquid  jet  by 
assuming  the  liquid  jet  and  the  pool  to  be  inviscid  liquids  and  the  gas  to  be  viscous. 
They  did  not  analyze  the  surface  depression. 

Bonetto,  et  al  [1993]  analyzed  an  inclined  plunging  liquid  jet.  They  assumed 
that  both  the  gas  and  the  liquid  were  inviscid  and  they  obtained  the  entrained  gas 
flow  rate.  The  volumetric  flow  rate  of  entrained  gas  is  one  of  the  most  important 
quantities  which  one  wants  to  compute  in  problems  of  this  type.  The  key  parameter 
in  such  evaluations  is  the  gas  gap  thickness  (ie,  the  shape  of  the  induced  surface 
depression). 

The  objective  of  this  work  was  to  evaluate  the  induced  surface  depression 

caused  by  a  plane  plunging  liquid  jet.  We  have  assumed  that  both  fluids  are 

inviscid  and  irrotational,  and  that  the  gas  region  is  at  constant  pressure.  Hence  the 

appropriate  equation  for  this  problem  is  Laplace's  equation  for  the  liquid  and  the 

associated  interfacial  jump  condition.  The  problem  may  be  described  by  two 

1  ^2 

parameters,  the  Weber  number.  We  = - — - and  the  Bond  number,  B^  =  — — —  , 

where  is  the  half-width  of  the  plunging  liquid  jet,  is  the  velocity  of  the  liquid 
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jet,  is  the  density  of  the  liquid,  g  is  the  gravitational  acceleration,  and  o  is  the 
surface  tension. 

We  nave  used  a  nonsingular  perturbation  technique  to  solve  the  problem. 
That  is,  we  have  expanded  the  solutions  for  relatively  small  We  number,  and 
substitute  these  expressions  into  Laplace's  equation  and  the  interfacial  jump 
condition.  Hquating  terms  of  the  same  order  in  We,  one  obtains  a  recursive  system 
of  equations,  which  we  have  solved  numerically  up  to  third  order. 


ASYMPTOTIC  EXPANSION 

The  purpose  of  this  section  is  to  compute  the  position  of  the  interface, 
y  =11  (x)  (shown  schematically  in  Figure-1). 

For  the  assumption  of  irrotational,  inviscid  flow  the  governing  equations  for 
the  liquid  field  are, 

v2y  =  0  (1) 


where  i]/  is  the  stream  function.  The  two  velocity  comporients  are  then  given  as. 


av 


(2) 


V 


aij/ 


(3) 


where  u  andv  are  the  velocity  components  along  the  x  and  y  axis,  respectively.  If 

aijr 

we  prescribe  a  value  of  \|/,  or  its  derivative  normal  to  the  surface,  for  every  point 

on  the  boundary,  we  have  a  w'ell-posed  problem.  As  shown  in  Fig.-2,  the  plane  x  =  0 
is  a  symmetry  plane.  Then  the  transverse  lateral  velocity,  u  must  vanish  for  every 
y  atx  =  0.  That  is,  from  Eq.  (2) 
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u(x  =  0  ,y)  (^  =  0  ,y)  =  0  (4) 

on  the  centerline  of  the  plunging  liquid  jet.  We  know  that  the  axial  velocity  must 
be  at  points  where  the  jet  is  impacting.  From  Fig-1  and  Eq.  (3),  we  see  that; 

V  (0  <  X  <  x^,y  =  0)  =  -  (0  <  X  <  XQ,y  =  0)  =  (5) 

For  X  >  Xq,  the  free  surface  position,  fj(x),  must  be  coincident  with  a  streamline. 
Without  loss  of  generality,  we  make  the  free  surface  coincident  with  the  streamline 
vj/  =  0.  Then, 

V(x,fi(^))=0  (6) 


We  assume  that  the  pressure  in  the  gas  region  is  constant  and  equal  to  zero.  The 
liquid  pressure  right  under  the  surface  is  related  to  the  curvature  of  the  surface  by: 


1  + 


dn. 


7;^  =  -p,  (x,ti(x)) 


Bernoulli's  equation  for  the  liquid  pressure  gives 


1 


-  p^  t<,fi(x))  =  +  2  [u^C'<,fl(x))  +  v^(x,f|(x))]  +  gilfx) 


(7) 


(8) 


Far  from  the  jet,  we  impose  the  boundary  conditions: 

\|/  (x  oo,y)  =  0  (9a) 


\|/  (x,  y  — >  oo)  =  0 


(9b) 


We  may  make  Eqs.  (l)-(9)  nondimensional,  using  as  the  length  scale  and,  as 
the  velocity  scale,  respectively.  Thus  we  have; 


o 

II 

r-j 

L> 

(10) 

with  boundary  conditions. 

o 

II 

o' 

II 

(11a) 

V|/  (x,  y^oo)  =  0 

(11b) 

y  (x->°o,y)  =  0 

(11:) 

(x,  y=0)  =  + 1 

0  <  X  <  1 

(lid) 

y(x,Ti  (x))  =  0 

1  <x 

(lie) 

where. 

<x 

II 

X 

II 

,  T]  =Tl/Xo 

u  =  u/V^ 

V  =  v/V( 

and  y  =  X  V 

The  velocity  components  may 

now  be  computed  from  the  stream 

function  using 

Eqs.  (2)  and  (3)  as. 

u(x,y)  =  (x,y) 

(12) 

v(x,y)  =  - (x,y) 

(13) 

The  interfacial  jump  condition  becomes 
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d^ri/dx“ 


{ 


.dn 


3/2 


1 

2 


[u^(x,n(x))  +  V^(x,Tl(x))]  +  T1 


(14) 


J 


Notice  that  the  two  multidimensional  parameters  are  the  Weber  number, 

1  ^  P^g^o 

We  =  2  P^  Bond  number,  =  — - —  ,  which  appear  in  the  boundary 

condition,  Eq.  (14), 

In  order  to  solve  the  problem  analytically  we  would  have  to  obtain  a  solution 
of  Laplace's  equation,  Eq.  (10),  that  satisfies  the  boundary  conditions,  Eq.  (11),  where 
the  free  surface  depression,  q(x),  comes  from  Eq.  (14).  We  note  that  the  solution  of 
Eq.  (14)  is  linked  to  the  solution  of  through  u  and  v;  thus  the  problem  is 
nonlinear. 

We  shall  solve  the  problem  for  small  We  using  a  perturbation  analysis.  First, 
we  may  expand  all  the  dependent  variables  in  terms  of  We  obtaining: 

n 

Ti(x)  =  Z  We‘  qj(x)  +  0(We"^b  (15a) 

i=l 

n 

\(/(x,y)  =  Z  We'  i(/i(^'y)  +  0(We''''’^)  (15b) 

i=:0 


Next  we  substitute  Eq.  (15b)  into  Eq.  (10)  to  obtain, 

+  WeV^v|/^  +  We^V^vi/^  +  . . .  =  0  (16) 

Equation  (16)  holds  for  all  We,  and  thus  all  V}/.  must  independently  satisfy  Laplace 
equations: 

=  0 


(17a) 
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=  0 

(17b) 

=  0 

(17c) 

Using  the  same  reasoning,  it  is  easy  to  show  that  all  the  homogeneous  boundary 
conditions,  Eqs.  (11a)  -  (11c),  leaH  to: 


9Yj 

ay  (x=o.y)  -  0 

(i>0) 

(18) 

Yj(x,y-^«>)  =  0 

(i>0) 

(19) 

Y.(x^oo,y)  =  0 

(i>C) 

(20) 

Let  us  next  expand  the  boundary  condition  in  Eq.  (lid): 

_  9\1/, 

(x,y =0)  +  We  (x,y =0)  +  We^  (x,y=0)  +  . . .  =  +1  (21 ) 

Obviously  the  zeroth  order  term  must  be  equal  to  +1  and  all  other  terms  must  be 
zero.  Thus, 

9\|/ 

(x,y=0)  =  +  1  (22) 

and 

d\\f. 

(x,y=0)  =  0  i  >  1  (23) 

The  boundary  condition  in  Eq.  (lie),  and  the  interfacial  jump  condition,  Eq.  (14), 
require  special  attention.  Using  Eq.  (15b),  Eq.  (lie)  can  be  rewritten  as: 

0  =  \|/(x,q  (x))  =  (x,Ti(x))  +  We  \(t^(x,Ti(x))  +  WeV2(>^''n(x))  +  ••• 


(24) 
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Performing  an  expansion  of  Vj/.in  terms  of  We  in  the  neighborhood  of  ri  =  0,  we 
obtain; 


0\j/  1 3“V|/  - 

0  =  (x,0)  +  (x,0)  n(x)  +  2  (x,0)  Ti‘^(x)  +  . . . 


+  We 


1 2 

'^1  n(x)  +  2  (X/O)  n  tx)  + . . . 


+ ... 


(25) 


Substituting  Eqs.  (15a)  and  (12)  into  Eq.  (25),  and  rearranging  (ie,  collecting  terms  of 
order  We‘),  we  obtain; 


+  0(We'^)  ... 


(26) 


Thus,  the  boundary  conditions  for  x  >  1  are, 
\|/^(x,0)  =  0 


(27a) 


\|/^(x,0)  =  -  u^q^ 


\|/2(x,0)  =  - 


1  2 


\|/2(x,0)  = 


1 


9u_ 


:i2 

1  5 


“l^2  ^2^1  2  ay  ’ll  ^  dy  ’ll’l2  +  6  ay2  ’ll 


(27b) 

(27c) 


(27d) 
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Notice  that  vj/j  has  homogeneous  boundary  conditions  for  i  >  1,  except  for  x  >  1, 
where  its  value  is  given  by  Eqs.  (27).  This  set  of  equations  is  the  result  of  the 
asymptotic  expansion  of  Eq.  (lie).  In  order  to  obtain  a  closed  set  of  equations,  let  us 
focus  our  attention  on  Eq.  (14),  the  interfacial  jump  condition.  We  first  expand  the 
left  hand  side  of  Eq.  (14)  and  then  the  right  hand  side. 

The  left  hand  side  of  Eq.  (14)  can  be  expanded  as: 


d“ri/dx“ 


dx" 


I  Ti,(x)  We'  +  0(We"^^) 

1=0 


f  .dn 
'  "  (d-x  > 


v.3/2 

) 


'  1  + 

/ 

d 

dx 

\ 

I  n.(x)  We' +  OlWe'"^^) 


1=0 


3/: 


d\  d'n,  d“Tio 

— ^  +  — =r  We  +  — f  We^  +  We"* 
dx"  dx“  d.x" 


f  d-ri3 

2  dxH 


dx“ 


+  We*^ 


dx 


2  dx  dx  '  2  dx- 


2„  > 


d“Tii  dn^  3d^n2fdTi,V  d“ri^ 


dx“ 


+  O(We^) 


(28) 


Similarly,  the  right  hand  side  of  Eq.  (14)  can  be  expanded  as, 
We[u^(x,q(x))  +  v^(x,q(x))]  = 


We  1 


du. 


u^(x,0)  +  (x,0)  q(x)  +  2  (x,0)  q%x)  +  . . 


0U. 


1 


+We(u^(x,0)  +"0y'  (x,0)q(x)  +  2  “^^(x,0)  q%x)  +  ...) 


3v, 


1^% 


Vq(x,0)  +  (x,0)  q(x)  +  2  (x,0)  q  (x)  + 
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( 


+  We 


V 


Vi(x,0)  ^  (x,0)  ri(x)  +  2“^^  (x,0)  il^x)  +  ... 


(29) 


Using  Eq.  (15a)  we  obtain: 

We[u“(x,q(x))  +  v^(x,q(x))] 


=  We 


Ou  ^ 

u^(x,o)  +  -^(x,o)  (n^  +  Wen,  +  wvn2  ■*■••• 


+  W  e 


Ou , 


u.(x,0) 


+  (x,o) (n^^  +  Wen,  +  vve“n2 ••M 


v^^(x,0)  +  (x,0)  (no('<)  +  Wen,  +  VVe“n2  ■*■•■•) 


+. . . 


+  We 


(X/0)  +  (x,0)  (n^  +  We  n,  +  We-n2  +...)  +  .. 


(29b) 


Finallv, 


11  ( 

(u"  +  v")  We  =  (u"  +v‘')  We  +  We^  ° 


2u. 


Ov 


rOv, 


+2v 


JJ 


+  We^  1 


1 

av,  av^  law^  2  9 

2 

[ay 

[3yn.-ayn,*2  3^2n,-vJ 

+  2v 


+  2u 


Ov 


1  ^-v,  2  ^^^2 


d\- 

o  o 


Ou 


1 


3“u,  2  1  0‘^u 


0  3 


11 


+  2 


aui 


la^u^  2 


2  ay^  '‘^ 


Hi  +  Li, 


(30a) 


Comparing  the  expression  for  the  right  hand  side  of  Eq.  (28)  with  the  right  hand  side 
of  Eq.  (30a)  we  obtain: 


d~qj  -) 

dx 


d-q. 

— 

dx" 


2  u 


au 


ay 


rn,  +u, 


(dx 


+  2  V, 


+  B„q 


0  'li 


(31a) 


(31b) 


(31c) 


The  boundary  conditions  for  these  ordinary  differential  equations  are: 


q  (x->oo)  =  0 


ciq, 

"dx  ~  ^ 


(32a) 

(32b) 


Equations  (27)  and  (31)  are  the  expanded  forms  of  Eqs.  (lie)  and  (14),  respectively. 
Equations  (17),  (18),  (19),  (20),  (22)  and  (23)  complete  the  set  of  equations  to  be  solved. 

ORDER  ZERO  SOLUTION 

In  the  previous  section  we  have  established  an  infinite  set  of  coupled 
equations  for  qj(x).  In  this  section  we  will  solve  the  order  zero  equation  analytically 

(i.e.,  for  the  functions  q^^,  V|/^,  u^,  v^).  The  equation  for  q^  is  Eq.  (31a), 

R 

^  +  Bo 


with  the  boundary  conditions 
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T1^(x->oo)  =  0 

(33a) 

dx 

(33b) 

The  solution  is 

=  0  (34) 

The  differential  equation  for  \}/^  is  Eq.  (17a), 

V-v|/^  =  0 

and  the  corresponding  boundary  conditions,  Eqs.  (18),  (19),  (20),  (22)  and  (23),  are: 


y^(x->oo,y)  =  0 

(35a) 

(35b) 

O 

II 

o 

II 

X 

(35c) 

8\|/ 

(x,y-0)  =  +  1 

0  <  X  <  1 

(35d) 

ax  “ 

X  >  1 

(35e) 

It  is  convenient  to  use  complex  variables  to  solve  this  two-dimensional  problem. 
Let  O  be  a  nondimensional  complex  analytic  function  given  by, 

0(z)  =  0(x,y)  +  i  \|/(x,y)  (36) 

where,  z  =  x  +  i  y,  and  0,  \|/  are  the  real  and  imaginary  parts  of  0(z),  respectively. 
The  function  0(z)  is  the  complex  potential,  ^  is  the  velocity  potential  function  and 
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\|/  is  the  corresponding  steam  function.  The  velocity  components  of  u,  v  are  given 
by. 


d<^  d\\i 
~  dx  ~  dz 

_  ^  _  dxjf 
~  dy  ~  '  dx 


(37a) 


(37b) 


=  u  - 1  V 


VVe  recall  the  boundary  conditions  (35d&e), 


-1  0  <  X  <  1 


A  function  of  x  and  y  that  satisfies  these  boundary  conditions  is: 

5^0  U  ("x-n  /^x+iV 

This  expression  is  the  imaginary  part  of  the  following  analytic  function 

dO  1  r  T 

-^  =  -[log  (z-l)-log  (z+1)] 


We  may  use  Eq.  (38)  to  calculate  the  velocity  components 
Up  =  Re((t>)  =  -  (log  [y“  +  (x-1)^]  -  log  [y^  +  (x+1)^]) 

.  W  (x+1)  (x-D) 

v^  =  -Im  (0)  =  -|^atg-y--atg— J 
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Table-I  shows  the  velocity  of  the  zeroth  order  solution  as  well  as  the  derivatives 
needed  for  solving  Eqs.  (27)  and  (31).  We  now  have  the  solution  of  the  problem  to 
order  zero. 

In  order  to  solve  the  problem  to  order-n,  we  must: 

(1)  Compute  u^j  (x,y=0)  and  VQ(x,y=0)  using  Eqs.  (42)  and  (43). 

(2)  Solve  Ti|(x)  from  (31b),  (32a),  (32b)  with  t1q(x,0),  Vq(x,0)  from  Eqs.  (42)  and  (43). 

(3)  Compute  V|/|(x,0)  from  Eq.  (27b)  for  x  >  1  using  q^fx),  Uq(x,0),  Vq(x,0). 

(4)  Solve  ^  =  0  using  yj(x,0)  to  obtain  \j/j(x,y). 

(5)  Compute  Ui(x,0)  =  (x,0),  Vi(x,0)  =  -  (x,0) . 

d“Ti2 

(6)  Evaluate  '  ,  2  using  Eq.  (31c). 

dx 

(7)  Compute  V|/2(x,0)  from  Eq.  (27c)  for  x  >  1,  using  q-,,  u^,  v^,  Up  Vp 

(8)  Solve  V^\y2  =  0  using  V|/2(X/0)- 

RESULTS 

Eigure-3  shows  the  analytical  results  for  the  zeroth  velocities  Vp,UQ  as  a 
function  of  the  lateral  position  x  for  different  axial  positions  y.  We  note  that,  as 
expected,  the  jet  is  spreading  as  we  move  downward.  Moreover,  the  singular  point 
at  X  =  1  in  the  lateral  velocity  is  clearly  visible  in  Eig-3(a). 

The  problem  consideed  so  far  is  spatially  unbounded  (x-  >  «>,y-  >  ««). 
However,  in  order  to  perform  the  numerical  integration  of  Eqs.  (17)  and  Eqs.  (31)  we 
need  finite  boundaries.  Thus,  we  have  taken  the  integration  domain  for  0  <  x  <  D,  0 
<  y  <  D.  In  the  computations  D  was  chosen  large  enough  to  have  the  velocities  at 
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the  boundaries  reduced  to  1%  of  the  maximum  value  V^.  The  corresponding 
physical  problem  is  a  jet  of  width  h  impacting  a  pool  of  width  2D  of  the  same  liquid. 
For  convenience  the  jet  contact  angle  was  7i/2,  and  Ti  is  measured  from  the  contact 
line  position.  That  is, 

ti(x=D)  =  0  (44a) 

^  ((x=D)  =  0  (44b) 


What  follows  is  a  detailed  description  of  Steps  (1)  -  (8)  from  the  previous 
section. 

Step  (1):  From  Eqs.  (42)  -  (43)  UQ(x,y=0)  must  be, 

u„(x,y=0)  =  ^  log  (45) 

and,  VQ(x,y=0)  =  0. 


Step  (2):  Equation  (31b)  results  in  the  following  ordinary  differential  equation, 


nv  -  B„iii  =  -  (u^  +  vj)  =  -  log-  ^ 


(46) 


for  the  interval  1  <  x  <  D  with  the  boundary  conditions. 


ri^(x=D)  =  0 

(47a) 

drii 

^'(x=D)  =  0 

(47b) 

Equation  (46)  has  an  integrable  singular  point  at  which  x=l.  We  have  numerically 
evaluated  Eq,  (46)  using  a  Runge-Kutta  (RK)  algorithm  with  adaptive  stepsize. 
Eigure-4(a)  shows  rij  as  a  function  of  x. 
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Since  Eq.  (46)  has  a  singular  point  at  x=l  we  have  performed  an  asymptotic 
analysis  of  x  near  x=l  to  verify  the  numerical  result.  Let  us  consider  the  solution  of 
Eq.  (46)  in  the  neighborhood  of  x=l.  An  asymptotic  expansion  of  the  right  hand  side 
of  Eq.  (46)  yields, 


1  ,  (x-1)  1  ,  ,  (x-1) 

Vlog2^~-^log2-^ 


(48) 


The  solution  which  is  valid  asymptotically  close  to  x=l  is. 


ri|(x)  =  q''^^  +  C2  e 


(49) 


The  main  contribution  for  the  integrals  on  the  right  hand  side  of  Eq.  (49)  comes 
from  small  values  of  (t-1).  Thus  w'e  can  expand  the  exponentials  as. 


With  this  analysis  we  have  shown  that  the  singularity  at  x=l  is  an  integrable 
one.  The  results  of  the  asymptotic  expansion  are  shown  in  Fig.  4  as  a  line  which 
essentially  coincides  with  the  numerical  result.  Indeed,  the  difference  between  the 
RK  solver  and  the  asymptotic  expansion  result  is  smaller  than  . 


\l/^(x,y==0)  is  computed  using  Eq.  (27b). 


\|/^(x,y=0)  =  -  Up(x,y=0)  ri^(x) 


(51) 
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The  resulting  boundary  condition  for  is  shown  in  Fig.  5(a).  Notice  that  has  a 
singularity  in  x=l  due  to  u^. 

Step  (4):  The  Laplace  equation  for  has  to  be  solved  with  the  boundary  conditions 
(BC)  given  in  Eqs.  (18)-(23).  We  numerically  solved  the  Laplace  equation  using  a 

five  points  finite  difference  scheme.  A  nine  points  scheme  was  also  implemented 
into  no  noticeable  improvement  over  the  five  point  scheme.  The  velocities  u^  v^ 

were  computed  using  a  differentiating  centered  first  order.  The  region  near  x  =  1 

(X“l) 

was  highly  refined  because  the  boundary  condition  y^  is  singular  there  (~  log  ~^)- 
We  have  verified  the  resulting  y^  using  an  asymptotic  expansion  of  Eq.  (52): 

y^(x,0)  =  -  u^q^  =  -  - - log(x-l)  =  y^^^’^(x,0)  ,  (1  <  x  <  D) 

y j  (x,0)  =  0  =  yf  ^""(x^O)  ,  (0  <  x  <  1 )  (53) 


An  analytic  continuation  of  Eq.  (53)  gives. 


asvm  Tl|(x"l)  3  i 

yi  ^  (z)  =  -  — ^  log  (z-1)  -  ^  log  [i(z-l)]  log(z-l) 


(54) 


The  real  part  of  Eq.  (54)  is  a  solution  of  the  Laplace  equation  and  has  the  same 
asymptotic  behavior  at  y.  near  x=l. 

Notice  that  although  y^(x,0)=  (x,0),  y^  becomes  significantly  different  than 

yasym  move  away  from  the  singularity,  the  results  near  x=l  are  almost 

coincident.  We  have  taken  advantage  of  the  linearity  of  the  Laplace  equation 
solving  numerically  for  the  difference  between  y^  and  That  is,  we  define 

diff,  ,  ,  .  asvm,  .  ,  ^ 

y^  (x,y)  =  yj(x,y)-y^  (x,y)  (55) 
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diff 


Notice  that  (x=l,y=0)  =  0  and  therefore  the  singularity  at  x=l  has  been  removed. 


diff 


The  boundary  conditions  for  are: 


\j/^‘^^(x,0)  =  \ir^(x,0)  -  \|/^^^(x,0) 

(56a) 

diff 

avi 

ax 

(56b) 

vfVy=0)  =  -vr’'"'<x.D) 

(56c) 

\j/j  (x=D,y)  =  -  (D,y) 

(56d) 

diff,  „  . 

(x=0,y)  =  0 

(56e) 

\|/^(x,y)  was  numerically  evaluated  using  the 

Laplacian  solver  with  Eqs.  (56).  We 

note  that  these  BCs  have  no  singularities. 

velocities  are  computed  using 

Finally,  the  stream  function  and  the 

'|/|(x,y)  =  v|;f‘^^(x,y)  +  \|/^^^’^(x,y) 

(57a) 

diff  asym 

u^(x,y)  =  U|  (x,y)  +  u^  (x,y) 

(57b) 

,  ,  diff,  .  asvm,  . 

v^(x,y)  =  v^  (x,y)  +  v^  ^  (x,y) 

(57c) 

The  maximum  difference  between  u^  computed  usign  Eqs.  (56)  and  (57),  and 

by  the  Laplacian  solver  applied  to  Eq.  (52),  was  less  than  2%  for  the  stream  function 
and  less  than  5%  for  the  velocities. 


Step  (5):  u^  and  v^  were  computed  from  the  stream  function  \j/^.  Figure  6  and  7 
show  both  velocity  iterates  as  a  function  of  the  position. 
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Step  (6):  We  calculated  ri2  from  Eq.  (31c)  using  the  RK  solver.  We  verified  the 

solution  in  the  neighborhood  of  x=l  as  described  previously.  The  result  is  shown  in 
Fig.  4.  We  see  that  the  x  interval,  where  is  significantly  different  than  zero,  is 
smaller  than  the  range  where  was  significantly  different  than  zero.  This  is  a 
general  result  of  our  calculations.  The  higher  the  order  of  the  iterate,  the  smaller 
the  range  where  the  values  of  the  iterate  are  significantly  different  from  zero. 

Step  (7):  We  compute  \if2(X/0)  using  Eq.  (31c).  The  result  is  shown  in  Fig.  5. 

Step  (8):  We  numerically  solved  for  \|r2(x,0)  using  the  Laplacian  solver.  We  verified 
the  numerical  result  as  described  for  Xj/j.  Notice  that  we  had  to  keep  track  of  three 
singularities  in  Eq.  (27c).  That  is. 


f  1  3^0  2 

V2(x,0)=-  Ujqj  +  Xoq2  + 

V  ^  J 


OU 

~  -  q|(x=l)  u^^^'^(x,0)  -  q2(x=l)  u^  ^'^(x,0)  +  ^|(x=l)  — ^ —  (x,0) 
asym  asvm  asym 

=  V21  +  ¥22'  +  ¥23 


then  we  obtained. 


asym  ,  .  asym,  .  asym,  ,  asym,  . 
\|/2  (z)  =  Y22  fz)  +  Y22  +  ¥23  fz) 


The  stream  functions  V2l'¥92'¥23  'vere  obtained  with  the  help  of  a  symbolic 

manipulator  an  because  of  their  complexity  the  results  will  not  be  presented  herein. 
As  before,  Y2 defined  as, 

diff  asym  .  . 

=  ¥2  ■  ¥2 


The  velocity  iterates  U2  V2  are  shown  in  Figs.  6&7. 


Step  (9):  The  third  order  iterates  rjj,  V3  were  computed  following  a  similar 

procedure.  The  results  are  shown  in  Figs.  (4)-(7).  We  see  that  significant  values  of 
the  third  order  iterates  are  restricted  to  a  smaller  range  than  the  second  order 
iterates. 

Figure  8  shows  the  position  of  the  interface  t]  as  a  function  of  the  lateral 
position  X  for  values  of  the  Weber  number  in  the  range  [0.1  to  0.8].  We  now  analyze 
the  behavior  of  the  ri(x)  as  we  move  from  the  rightmost  position  in  the  plot  (x=l.l) 
to  lower  values  of  x.  For  We  =  0.01,  T|(x)  increases  as  x  approaches  the  jet  impact 
point  X  =  1  without  any  local  minimum  (ie,  no  surface  depression  is  preseid).  For 

We  =  0.8,  on  the  other  hand,  as  x  is  decreased,  we  first  have  a  local  minimum  in  T], 
(T^maX),  ^  local  minimum  depression  is  present).  We  define  the 

depth  of  the  depression,  d,  as: 


,  max  min 

d  =  n  -n 


(61) 


Figure  9  shows  the  depth  of  the  depression  d  as  a  function  of  the  Weber  number. 
There  is  no  depression  (d=0)  up  to  a  critical  Weber  number  (We^.).  For  We  larger 

than  the  critical  value,  the  depth  of  the  depression  increases  rapidly  with  the  We. 

The  depression  width  6  is  defined  as  the  gas  gap  corresponding  to  a 
depression  of  d/2.  Figure  10  shows  the  depression  width  5  as  a  function  of  the  We 
number.  We  see  that  for  high  We  numbers  5  levels  off  to  a  value  of  about  0.013. 

Finally,  we  present  a  set  of  results  where  the  effect  of  gravity  has  been 
investigated.  In  Figs.  11-17  the  Bond  number  is  equal  to  Bo  =  1.  The  interface 
position  iterates,  Tii(Bo  =  1),  have  the  same  qualitative  shape  compared  to  Tij(Bo  =  0). 

However,  the  magnitude  of  the  iterates  is  larger  in  the  case  with  Bo  =  1.  The  same  is 
true  for  the  stream  function  and  velocities  iterates.  This  effect  is  a  direct 
consequence  of  the  larger  values  of  the  qj.  Notice  that  u^  has  the  same  value  in  the 
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Bo  =  0  and  Bo=  1  calculations.  Thus,  if  Iti^(Bo  =  1)1  is  greater  than  It||(Bo  =  0)  I  the 
same  relationship  will  hold  between  ly^(Bo  =  1)  I  and  lx|/|(Bo  =  0)1. 

Figure  15  shows  the  position  of  the  interface  as  a  function  of  the  lateral 
position  for  different  We  number.  The  most  noticeable  effect  of  the  gravity  is  to 
flatten  the  interface. 

Figure  16  shows  the  depth  of  the  depression  (d)  as  a  function  of  the  We 
number  for  Bo  =  1.  The  interface  has  a  larger  depression  for  Bo  =  1  than  Bo  =  0. 
Figure  17  shows  the  width  of  the  depression  as  a  function  of  the  We  number  for  Bo 
=  1.  Surprisingly,  even  though  the  shape  of  the  interfaces  are  different  than  the  Bo  = 
0  the  width  is  nearly  equal. 

SUMMARY  AND  CONCLUSIONS 

This  paper  presents  the  results  of  an  analysis  based  on  a  nonsingular 
perturbative  technique  for  the  surface  depression  produced  by  a  plunging  liquid  jet. 
The  first  three  terms  of  a  Taylor  series  expansion  in  the  Weber  number  have  been 
obtained  for  the  surface  depression.  This  approximation  of  the  surface  depression 
gives  correct  values  for  small  and  moderate  Weber  numbers  (ie.  We  <  1).  Flowever, 
for  We  greater  than  unity,  higher  order  terms  become  important  and  the  analysis 
presented  here  is  no  longer  valid.  The  Bond  number  appears  as  a  parameter  in  the 
system  of  equations.  No  expansion  is  necessary  for  the  Bond  number.  The  analysis 
presented  here  is  valid  for  all  Bo  numbers. 

The  results  described  in  this  paper  show  that  as  the  Weber  number  is 
increased,  the  terms  that  gain  importance  (i.e.,  the  terms  of  higher  order)  correspond 
to  a  surface  depression  that  is  increasingly  narrower  in  the  horizontal  direction.  For 
a  Weber  number  of  the  order  of  unity  the  surface  tension  is  no  longer  strong 
enough  to  keep  the  system  stable  and  an  instability  leads  to  air  entrainment.  That  is, 
for  values  of  the  surface  tension  going  to  infinity  (i.e.,  Weber  number  going  to  zero). 
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the  slope  of  the  surface  is  very  small,  however,  as  the  Weber  number  is  increased 
the  slope  also  increases.  For  a  critical  value  of  the  Weber  number,  the  slope  of  the 
surface  is  such  that  the  surface  tension  is  not  large  enough  to  keep  the  pool  surface 
from  getting  near  the  plunging  liquid  jet  and  thus  air  entrainment  is  produced.  We 
have  assumed  in  this  paper  that  the  position  of  the  interface  can  be  written  as  r|(x), 
that  is  a  steady  state  exists  where  the  position  of  the  interface  is  a  single-valued 
function  of  the  lateral  position.  For  the  critical  value  of  the  Weber  number  that 
corresponds  to  the  threshold  for  air  entrainment  the  interface  is  not  a  single-valued 
and  also  it  is  likely  the  air  entrainment  process  involves  a  transient  phenomena. 
Thus  we  believe  that  our  analysis  is  not  appropriate  to  compute  the  threshold  value 
of  the  Weber  number  for  air  entrainment.  However,  we  are  able  to  describe  the 
route  to  air  entrainment  qualitatevily.  We  have  seen  that  the  higher  the  order  of 
the  position  of  the  interface  iterate,  the  larger  the  maximum  slope  of  the  iterate.  For 
small  We  numbers  the  low  order  iterates  dominate.  As  the  We  number  is  increased, 
the  slope  of  the  interface  is  increased  because  higher  order  terms  dominate.  For  a 
value  of  We  the  interface  has  a  very  large  slope  (tangent  of  the  interface  almost 
vertical).  For  a  We>Wet  the  inertia  forces  overpower  the  capillary  forces  and  the 
interface  is  not  stable  any  more. 

The  results  presented  in  this  paper  cannot  predict  the  route  to  instability  that 
produces  air  entrainment  for  two  reasons.  First,  it  is  not  clear  that  the 
approximation  of  the  surface  depression  is  valid  for  Weber  numbers  that  are  high 
enough  to  produce  air  entrainment.  Second,  the  shape  of  the  surface,  ti(x),  has  to  be 
a  monovalued  function  of  x  and  it  appears  that  at  the  point  where  air  entrainment 
is  produced,  for  every  x  (horizontal  position)  there  is  more  than  one  ri  (vertical 
position  of  the  surface  depression). 
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It  appears  that  it  would  be  useful  to  compute  the  surface  position  using  an 
appropriate  multidimensional  Computational  Fluid  Dynamics  (CFD)  tool  having 
surface  tracking  capability.  This  will  allow  the  relaxation  of  the  small  and  moderate 
Weber  numbers  assumption. 
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TABLE-I 


SUMMARY  OF  THE  ORDER-ZERO  SOLUTIONS 


1  x-1 

Uo(x.0)=-log^ 


Vq  (x,0)  =  0 


(x,0)-0 


^y  7r(x-l)(x+l) 

4x 

dy  n(x-l)  (x+1) 


T  (x,0)  =  0 


— f  (x,0)  =  0 

dV' 


^  4(3x^-H) 

"  ■  rtx-l)^(x+I)3 


1  ry-  +  (x-l)2 

%(x.y)  =  2„Iog|^y2V,„i,2 

1  (  X  +  1  X-1 
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FIGURE  CAPTIONS 

Figure  1  Illustration  of  the  surface  depression  produced  by  a  plunging 
liquid  jet. 

Figure  2  Boundary  conditions  for  the  stream  function 
y  =  oo  Vj/  =  0  \j/  =  0 

X  =  oo  y  X  Xq  Vj 

symmetry  plane 

Figure  3(a)  Solution  of  the  transverse  velocity  of 

order  zero  for  different  axial  positions:  (a)  y=0  (undisturbed 
pool  surface  level),  (b)  y=0.5,  (c)y=l.,  and  (d)  y=2. 

(a)  (b)  (c)  (d)  X 

Figure  3(b)  Solution  of  the  axial  velocity  of  order 

zero  for  different  axial  positions:  (a)  y=0  (undisturbed  pool 
surface  level),  (b)  y=0.5,  (c)y=l.,  and  (d)  y=2. 

(a)  (b)  (c)  (d)  X 

Figure  4  Iterates  of  the  surface  position,  rijfx),  corresponding  to  Bo=0 

’ll 
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^3 

X 

Figure  5  Iterates  of  the  surface  stream  function,  \j/j(x,y=0), 
corresponding  to  Bo=0 

^2 

^3 

X 

Figure  6  Iterates  of  the  lateral  velocity,  Ui(x,y=0),  corresponding  to 
Bo=0 


^2 

X 


Iterates  of  the  axial  velocity,  Vj(x,y=0),  corresponding  to  Bo=0 


''2 


Figure  7 
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X 

Figure  8  Shape  of  the  interface  for  different  Weber  numbers(Bo=0) 

Tj  X  We  =  0.1  We  =  0.8 

Figure  9  Gas  gap  depth  d  as  a  function  of  the  Weber  number  (Bo=0). 

Figure  13  Definition  of  the  gas  gap  width,  6 

Liquid 

Jet 

Figure  10  Gas  gap  width  5  as  a  function  of  the  Weber  number  (Bo=0) 

5  We 

Figure  1 1  Iterates  of  the  surface  position,  Tjjfx),  corresponding  to  Bo=l 

’ll 

’12 

’I3 

X 

Figure  12  Iterates  of  the  surface  stream  function,  vj/j(x,y=0), 
corresponding  to  Bo=l 


Figure  13  Iterates  of  the  lateral  velocity,  Uj(x,y=0),  corresponding  to 
Bo=l 


Figure  14  Iterates  of  the  axial  velocity,  Vi(x,y=0),  corresponding  to  Bo=l 


Figure  15  Shape  of  the  interface  for  different  Weber  numbers(Bo=l) 
71  X  We  =  0.1  We  =  0.8 
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Figure  16  Gas  gap  depth  d  as  a  function  of  the  Weber  number  (Bo=l). 
Figure  17  Gas  gap  width  6  as  a  function  of  the  Weber  number  (Bo=l) 
5  We 


Figure  1  Illustration  of  the  surface  depression  produced  by  a  plunging 
liquid  jei. 


symm 


undary  conditions  for  the  stream  function 


Solution  of  the  transverse  velocity  of 
order  zero  for  different  axial  positions;  (a)  y=0  (undisturbed 
pool  surface  level),  (b)  y=0.5,  (c)y=l.,  and  (d)  y=2. 


Figuro  3(b)  Solution  of  the  axial  velocity  of  order 

zero  for  different  axial  positions:  (a)  y=0  (undisturbed  pool 
surface  level),  (b)  y=0.5,  (c)y=l.,  and  (d)  y-l. 


Figure  5  Iterates  of  the  surface  stream  function,  \j/j(x,y=0), 
corresponding  tc  Bo=0 


Figure  6  Iterates  of  the  lateral  velocity,  Ui(x,y=0),  corresponding  to 


o 


o 


Figure  7  Iterates  of  the  axial  velocity,  Vj(x,y=0),  corresponding  to  Bo- 


Figure  12  Iterates  of  the  surface  stream  function,  yi(x,y=0), 
corresponding  to  Bo=l 


Figure  14  Iterates  of  the  axial  velocity,  vj(x,y=0),  corresponding  to 


